package scu.maqiang.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q41;
import scu.maqiang.fes.FES2Q81;
import scu.maqiang.mesh.Mesh2Q4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DoubleIntegerPair;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;
import static scu.maqiang.numeric.Constants.*;

import java.util.Arrays;

public class SOTS2DQ481 {
	public SOTS2DQ481(FES2Q41 fsQ4, FES2Q81 fsQ8) {
		this.fsQ4 = fsQ4;
		this.fsQ8 = fsQ8;
	}
	
	/**
	 * 求解二阶热传导系数的一阶单胞问题，均匀化解以及二阶单胞问题
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 * @param SOCS
	 */
	public void KxxCellSolution(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D, double[][] SOCS) {
		SRMatrix A = new SRMatrix(fsQ4.GetNdof());
		fsQ4.assembleStiff(kxxFunc, param, BVPType.COMMON, A);
		double[][] RHS = new double[2][fsQ4.GetNdof()];
		fsQ4.assembleGlobalVector(kxxFunc, param, this::elementkxxFirstOrder, BVPType.COMMON, RHS);
		fsQ4.applyBC_MBN(A, Direct.All, 1, 2, 3, 4);
		fsQ4.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4);
		fsQ4.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4);
		NewIterSSolver solver = new NewIterSSolver(A);
		MVO.fill(FOCS, 0.0);
		solver.PCGSSOR(RHS[0], FOCS[0], 1.5, 1);
		solver.PCGSSOR(RHS[1], FOCS[1], 1.5, 1);
		//System.exit(0);
		
		computeHomogenizedKxx(kxxFunc, param, FOCS, D);
		
		RHS = new double[4][fsQ4.GetNdof()];
		fsQ4.assembleGlobalVector(kxxFunc, param, FOCS, D, this::elementKxxSecondOrder, BVPType.COMMON, RHS);
		MVO.fill(SOCS, 0.0);
		fsQ4.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4);
		fsQ4.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4);
		fsQ4.applyBC_MBN(RHS[2], Direct.All, 0.0, 1, 2, 3, 4);
		fsQ4.applyBC_MBN(RHS[3], Direct.All, 0.0, 1, 2, 3, 4);
		solver.PCGSSOR(RHS[0], SOCS[0], 1.5, 1);
		solver.PCGSSOR(RHS[1], SOCS[1], 1.5, 1);
		solver.PCGSSOR(RHS[2], SOCS[2], 1.5, 1);
		solver.PCGSSOR(RHS[3], SOCS[3], 1.5, 1);
	}
	
	/**
	 * �����ȴ���ϵ��һ�׵��������Ҷ���������
	 * @param xy
	 * @param coef
	 * @param tp
	 * @param eleXY
	 */
	public void elementkxxFirstOrder(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] eleXY) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] v = new double[4];
		double[] dx = new double[4];
		double[] dy = new double[4];
		Arrays.fill(eleXY[0], 0.0);
		Arrays.fill(eleXY[1], 0.0);
		for(int i = 0; i < 2; i++) {
			double xi = GP2[i];
			for(int j = 0; j < 2; j++) {
				double et = GP2[j];
				double detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
				double c = MVO.dot_product(coef[0], v);
				double wxarea = GW2[i] * GW2[j] * detJ * c;
				for(int k = 0; k < 4; k++) {
					eleXY[0][k] += -dx[k] * wxarea;
					eleXY[1][k] += -dy[k] * wxarea;
				}
			}
		}
	}
	
	public void elementKxxMatrixFirstOrder(int[] ele, double[][] xy, double[][] mat, BVPType tp, double[][] eleXY) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] v = new double[4];
		double[] dx = new double[4];
		double[] dy = new double[4];
		Arrays.fill(eleXY[0], 0.0);
		Arrays.fill(eleXY[1], 0.0);
		for(int i = 0; i < 2; i++) {
			double xi = GP2[i];
			for(int j = 0; j < 2; j++) {
				double et = GP2[j];
				double detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
				double wxarea = GW2[i] * GW2[j] * detJ;
				for(int k = 0; k < 4; k++) {
					eleXY[0][k] += -(mat[0][0] * dx[k] + mat[1][0] * dy[k]) * wxarea;
					eleXY[1][k] += -(mat[0][1] * dx[k] + mat[1][1] * dy[k]) * wxarea;
				}
			}
		}
	}
	
	/**
	 * ��һ�׵������������ȴ������Ȼ�ϵ������
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 */
	public void computeHomogenizedKxx(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D) {
		double[][] xy = new double[4][2];
		Mesh2Q4 mesh = (Mesh2Q4)fsQ4.getMesh();
		double[] NN1 = new double[4];
		double[] NN2 = new double[4];
		double wholeArea = 0.0;
		double[][] coef = new double[kxxFunc.length][4];
		MVO.fill(D, 0.0);
		double[] dx = new double[4];
		double[] dy = new double[4];
		double[] v = new double[4];
		double xi, et;
		for(int k = 0; k < mesh.getNt(); k++) {
			int[] ele = mesh.getElement(k);
			mesh.getCoordPerElement(ele, xy);
			MVO.getSlice(FOCS[0], ele, NN1);
			MVO.getSlice(FOCS[1], ele, NN2);
			double[] x = MVO.col(xy, 0);
			double[] y = MVO.col(xy, 1);
			fsQ4.getCoefFromFunc(kxxFunc, xy, mesh.getElementLabel(k), param, coef);
			for(int i = 0; i < 2; i++) {
				xi = GP2[i];
				for(int j = 0; j < 2; j++) {
					et = GP2[j];
					double detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
					double mKxx = MVO.dot_product(coef[0], v);
					double wxarea = GW2[i] * GW2[j] * detJ * mKxx;
					D[0][0] += (1 + MVO.dot_product(NN1, dx)) * wxarea;
					D[0][1] += MVO.dot_product(NN2, dx) * wxarea;
					D[1][0] += MVO.dot_product(NN1, dy) * wxarea;
					D[1][1] += (1 + MVO.dot_product(NN2, dy)) * wxarea;
				}
			}
			wholeArea += CG.area2DQ4(xy);
		}
		MVO.mul(D, 1.0 / wholeArea);
	}
	
	/**
	 * �����ȴ���ϵ�����׵����������Ҷ���������
	 * @param xy
	 * @param kxx
	 * @param FOCS
	 * @param homo
	 * @param tp
	 * @param eleXY
	 */
	public void elementKxxSecondOrder(int[] ele, double[][] xy, double[][] kxx, double[][] FOCS, double[][] homo, BVPType tp, double[][] eleXY) {
		
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] v = new double[4];
		double[] dx = new double[4];
		double[] dy = new double[4];
		Arrays.fill(eleXY[0], 0.0);
		Arrays.fill(eleXY[1], 0.0);
		Arrays.fill(eleXY[2], 0.0);
		Arrays.fill(eleXY[3], 0.0);
		for(int i = 0; i < 2; i++) {
			double xi = GP2[i];
			for(int j = 0; j < 2; j++) {
				double et = GP2[j];
				double detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
				double mkxx = MVO.dot_product(kxx[0], v);
				double wxarea = GW2[i] * GW2[j] * detJ;
				for(int k = 0; k < 4; k++) {
					eleXY[0][k] += ((mkxx - homo[0][0] + mkxx * MVO.dot_product(dx, FOCS[0])) * v[k] - mkxx * MVO.dot_product(v, FOCS[0]) * dx[k]) * wxarea;
					eleXY[1][k] += ((     - homo[0][1] + mkxx * MVO.dot_product(dy, FOCS[0])) * v[k] - mkxx * MVO.dot_product(v, FOCS[0]) * dy[k]) * wxarea;
					eleXY[2][k] += ((     - homo[1][0] + mkxx * MVO.dot_product(dx, FOCS[1])) * v[k] - mkxx * MVO.dot_product(v, FOCS[1]) * dx[k]) * wxarea;
					eleXY[3][k] += ((mkxx - homo[1][1] + mkxx * MVO.dot_product(dy, FOCS[1])) * v[k] - mkxx * MVO.dot_product(v, FOCS[1]) * dy[k]) * wxarea;
				}
			}
		}
	}
	

	
	public void assembleHeatSOTSSolution_Square(FES2Q41 fsFine, double[] T_0, double[] T_1, double[] T_2,
                                                FES2Q81 fsHomo, double[] T_Homo, FES2Q41 fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		double[] xy = new double[2];
		double[] cellXY = new double[2];
		int[] eleIdx8 = new int[8];
		double[] eleVal8 = new double[8];
		int[] eleIdx4 = new int[4];
		double[] eleVal4 = new double[4];
		double[][] grad = new double[2][fsHomo.GetNdof()];
		double[][] Se = new double[3][fsHomo.GetNdof()];
//		fsHomo.computeGradient2(T_Homo, grad);
//		fsHomo.computeSecondDerivative(T_Homo, Se);
		fsHomo.computeFirstAndSecondDerivative(T_Homo, grad, Se);
		double val = 0;
		double corrector1, corrector2;
		
		double[][] homoCoord = new double[8][2];
		double[][] homoCoef = new double[8][8];
		fsHomo.mesh.getCoordPerElement(fsHomo.mesh.getElement(0), homoCoord);
		CG.Q8ShapeFunctionCoef(homoCoord, homoCoef);
		DoubleIntegerPair homoData = new DoubleIntegerPair(0.0, 0);

		double[][] cellCoord = new double[4][2];
		double[][] cellCoef = new double[4][4];
		fsCell.mesh.getCoordPerElement(fsCell.mesh.getElement(0), cellCoord);
		CG.Q4ShapeFuncionCoef(cellCoord, cellCoef);
		DoubleIntegerPair cellData = new DoubleIntegerPair(0.0, 0);
		
		for(int i = 0, n = fsFine.GetNdof(); i < n; i++) {
			fsFine.mesh.getCoordPerNode(i, xy);
			cellXY[0] = xy[0] * N - (int)(xy[0] * N);
			cellXY[1] = xy[1] * N - (int)(xy[1] * N);
			
			//val = fsHomo.getValue(T_Homo, xy, eleIdx8, eleVal8);
			//T_0[i] = val;
			
//			Pair<Double, Integer> homodata = fsHomo.getValue(T_Homo, xy, previousHomoEleIdx, homoCoef, eleIdx8, eleVal8);
//			T_0[i] = homodata.getFirst();
//			previousHomoEleIdx = homodata.getSecond();
			fsHomo.getValue_Square(T_Homo, xy, homoData, homoCoef, eleIdx8, eleVal8);
			T_0[i] = homoData.first;
			//System.out.println("xy:" + Arrays.toString(xy) + "\t" + Arrays.toString(eleIdx8) + "---" + val + "---" + Arrays.toString(eleVal8));
			//try {
			//	Thread.sleep(200);
			//} catch (InterruptedException e) {
				// TODO Auto-generated catch block
			//	e.printStackTrace();
			//}
			
			fsCell.getValue_Square(FOCS[0], cellXY, cellData, cellCoef, eleIdx4, eleVal4);
			val = cellData.first;
			
			//Pair<Double, Integer> celldata = fsCell.getValue(FOCS[0], cellXY, previousCellEleIdx, cellCoef, eleIdx4, eleVal4);
			//val = celldata.getFirst();
			//previousCellEleIdx = celldata.second;
			
			//val = fsCell.getValue(FOCS[0], cellXY, eleIdx4, eleVal4);
			corrector1 = 0.0;
			
			corrector1 += val * MVO.dot_product(grad[0], eleIdx8, eleVal8);
			corrector1 += MVO.dot_product(FOCS[1], eleIdx4, eleVal4) * MVO.dot_product(grad[1], eleIdx8, eleVal8);
			
			T_1[i] = T_0[i] + ep * corrector1;
			corrector2 = 0.0;
			
			corrector2 +=  MVO.dot_product(SOCS[0], eleIdx4, eleVal4) * MVO.dot_product(Se[0], eleIdx8, eleVal8);
			corrector2 +=  MVO.dot_product(SOCS[1], eleIdx4, eleVal4) * MVO.dot_product(Se[1], eleIdx8, eleVal8);
			corrector2 +=  MVO.dot_product(SOCS[2], eleIdx4, eleVal4) * MVO.dot_product(Se[1], eleIdx8, eleVal8);
			corrector2 +=  MVO.dot_product(SOCS[3], eleIdx4, eleVal4) * MVO.dot_product(Se[2], eleIdx8, eleVal8);
			T_2[i] = T_1[i] + ep  * ep * corrector2;
		}
	}

	FES2Q41 fsQ4;
	FES2Q81 fsQ8;
}
